function [k_c, Ma_c] = neutral_stab_curves_naive

N = 20;
N_s = 20;
Ma_c = zeros(N, 1);
k_c = zeros(N, 1);


k_hat = zeros(N_s, 1);
ev_hat = k_hat;

p = params;


Ma = logspace(6, 4, N);
Ma_c = Ma;


for n = 1:N
    
    p.Ma = Ma(n);
    if n == 1
        k_hat(1) = 1e-1;
    else
        k_hat(1) = 1e-1;
    end
    
    ev_hat(1) = comp_eigs(k_hat(1), p);
    
    k_hat(2) = k_hat(1) + 1e-5;
    ev_hat(2) = comp_eigs(k_hat(2), p);
    
    d_ev = (ev_hat(2) - ev_hat(1)) / (k_hat(2) - k_hat(1));
    
    for its = 3:N_s
        
        
        k_hat(its) = k_hat(its-1) - ev_hat(its-1) / d_ev;
        ev_hat(its) = comp_eigs(k_hat(its), p);
        
        fprintf('%.4e\n', abs(ev_hat(its)));
        
        if (abs(ev_hat(its)) < 1e-10)
            k_c(n) = k_hat(its);
            fprintf('---------------------------\n');
            break
        end
        
        d_ev = (ev_hat(its) - ev_hat(its-1)) / (k_hat(its) - k_hat(its-1));
        
    end
    
    if (its == N_s )
        fprintf('did not converge!\n');
        k_c = k_c(1:n-1);
        Ma_c = Ma_c(1:n-1);
        break;
    end
    
end


loglog(k_c, Ma_c, 'k', k_c, 80 * (1 / p.delta + (1 - 2 * p.beta) ./ k_c.^2) / p.beta / (1 - p.beta),'r*');